//
// Created by Yurzi on 2021/3/21.
//

#include "RegressionPrediction.h"
#include <math.h>

//回归参数
double a[DEGREE];

//判断2浮点数大小
int doublecmp(double x,double y){
    if(fabs(x-y)<=MISS){
        return 0;
    } else if(x>y){
        return 1;
    } else{
        return -1;
    }
}

//生成以原料id为索引的统计数据的链表
BinaryPair* initListFromMTById(LogNode *rear,int MtId){
    BinaryPair *Bhead=NULL;
    BinaryPair *Brear=NULL;
    BinaryPair *Bp=NULL;
    BinaryPair *Bp0=NULL;
    LogNode *left=rear;
    LogNode *right=rear;
    LogNode *Lp=rear;
    int count=0;
    //遍历链表
    while (left!=NULL&&right!=NULL&&count<DEGREE*3){
        //维护一日区间
        while ((left!=NULL)&&(((right->time)/3600/24)-((left->time)/3600/24))<=0){
            left=left->prv;
        }
        Bp=(BinaryPair*) malloc(sizeof(BinaryPair));
        Bp->prv=NULL;
        Bp->next=NULL;
        Bp->x=count;
        Bp->y=0;
        if(count==0){
            Bhead=Bp;
            Brear=Bp;
        }
        if(count>0){
            Bp->prv=Bp0;
            Bp0->next=Bp;
        }
        Bp0=Bp;
        ++count;//计数加一
        Lp=right;
        while (Lp!=left){
            for (int i = 0; i <MAX_CARSIZE ; ++i) {
                if(Lp->content.materials[i].materialId==MtId){
                    Bp->y+=Lp->content.materials[i].amount;
                }
            }
            Lp=Lp->prv;
        }

        right=left;
    }
    Brear=Bp;
    return Brear;
}

//逆序
BinaryPair *reverseBList(BinaryPair *rear){
    BinaryPair *head=rear;
    BinaryPair *p=rear;
    BinaryPair *p0=NULL;
    int flag=0;
    //双链表逆序
    while (p!=NULL){
        p->next=p->prv;
        p->prv=p0;
        p0=p;
        p=p->next;
        flag=1;
    }
    if (flag){
        p0=head;
        int count=1;
        while (p0!=NULL){
            p0->x=count++;
            p0=p0->next;
        }
    }
    return head;
}

//初始化矩阵
void initMatrix(double **X_matrix,double **Y_matrix,int line,BinaryPair *list){
    if(line==0)return;
    BinaryPair *ptr=list;
    for (int i = 0; i <line ; ++i) {
        if(ptr!=NULL) {
            Y_matrix[i][0] = ptr->y;
        }else{
            Y_matrix[i][0]=1;
        }
        for (int j=0; j < DEGREE; j++){
            if(ptr!=NULL) {
                X_matrix[i][j] = pow(ptr->x, j + 1);
            }else{
                X_matrix[i][j]=1;
            }
        }
        if(ptr!=NULL)ptr=ptr->next;
    }
}

//矩阵乘法
void matrix_M(double **res_matrix,double **a_matrix,double **b_matrix,int line1,int row2,int middle){
    int i=0;
    int j=0;

    double temp=0;
    for ( i= 0;  i<line1; ++i) {
        for (j= 0; j < row2; ++j) {
            for (int k=0;k<middle;++k){
                temp+=a_matrix[i][k]*b_matrix[k][j];
            }
            res_matrix[i][j]=temp;
            temp=0;
        }
    }
}

//矩阵转置
void matrix_T(double **res_matrix,double **matrix,int line,int row){
    //res转置后矩阵
    //matrix转置前矩阵
    //line 行数
    //row 列数
    int i; //行
    int j; //列
    double temp=0;
    for (i=0;i<line;i++){
        for (j=0;j<row;j++){
            res_matrix[j][i]=matrix[i][j];
        }
    }
}

//矩阵求逆
int  matrix_inv(double **matrix, int ndimen){
    //	matrix:矩阵
    // 	ndimen :维数
    double tmp, tmp2, b_tmp[20], c_tmp[20];
    int k, k1, k2, k3, j, i, j2, i2, kme[20], kmf[20];
    i2 = j2 = 0;

    for (k = 0; k < ndimen; k++)
    {
        tmp2 = 0.0;
        for (i = k; i < ndimen; i++)
        {
            for (j = k; j < ndimen; j++)
            {
                if (fabs(matrix[i][j] ) <= fabs(tmp2))
                    continue;
                tmp2 = matrix[i][j];
                i2 = i;
                j2 = j;
            }
        }
        if (i2 != k)
        {
            for (j = 0; j < ndimen; j++)
            {
                tmp = matrix[i2][j];
                matrix[i2][j] = matrix[k][j];
                matrix[k][j] = tmp;
            }
        }
        if (j2 != k)
        {
            for (i = 0; i < ndimen; i++)
            {
                tmp = matrix[i][j2];
                matrix[i][j2] = matrix[i][k];
                matrix[i][k] = tmp;
            }
        }
        kme[k] = i2;
        kmf[k] = j2;
        for (j = 0; j < ndimen; j++)
        {
            if (j == k)
            {
                b_tmp[j] = 1.0 / tmp2;
                c_tmp[j] = 1.0;
            }
            else
            {
                b_tmp[j] = -matrix[k][j] / tmp2;
                c_tmp[j] = matrix[j][k];
            }
            matrix[k][j] = 0.0;
            matrix[j][k] = 0.0;
        }
        for (i = 0; i < ndimen; i++)
        {
            for (j = 0; j < ndimen; j++)
            {
                matrix[i][j] = matrix[i][j] + c_tmp[i] * b_tmp[j];
            }
        }
    }
    for (k3 = 0; k3 < ndimen;  k3++)
    {
        k  = ndimen - k3 - 1;
        k1 = kme[k];
        k2 = kmf[k];
        if (k1 != k)
        {
            for (i = 0; i < ndimen; i++)
            {
                tmp = matrix[i][k1];
                matrix[i][k1] = matrix[i][k];
                matrix[i][k] = tmp;
            }
        }
        if (k2 != k)
        {
            for(j = 0; j < ndimen; j++)
            {
                tmp = matrix[k2][j];
                matrix[k2][j] = matrix[k][j];
                matrix[k][j] = tmp;
            }
        }
    }
    return (0);
}

//拟合函数
void fitingFunc(BinaryPair *list){
    if(list==0)return;
    //构建矩阵
    //得到矩阵的行数
    int line=logDate_Length(list);
    //定义X矩阵与Y矩阵
    double **X_matrix=init2DMatrix(line, DEGREE);
    double **Y_matrix=init2DMatrix(line,1);
    double **res_matrix=init2DMatrix(DEGREE, 1);

    //输入矩阵内容
    initMatrix(X_matrix,Y_matrix,line,list);

    //求最小二乘解 使用最小二乘解求出最优拟合参数
    //求X的转置
    double **X_matrix_T=init2DMatrix(DEGREE, line);
    matrix_T(X_matrix_T, X_matrix, line, DEGREE);

    double **Temp1_matrix=init2DMatrix(DEGREE, DEGREE);
    matrix_M(Temp1_matrix, X_matrix_T, X_matrix, DEGREE, DEGREE, line);


    matrix_inv(Temp1_matrix, DEGREE); //因为X的值为天数不可能为0，变相的矩阵定是满秩矩阵

    double **Temp2_matrix=init2DMatrix(DEGREE, line);
    matrix_M(Temp2_matrix, Temp1_matrix, X_matrix_T, DEGREE, line, DEGREE);

    matrix_M(res_matrix, Temp2_matrix, Y_matrix,DEGREE,1, line);

    //将结果保存
    for(int i=0; i < DEGREE; i++){
        a[i]=(res_matrix[i][0]);
    }

    //释放空间
    free2DMatrix(X_matrix,line);
    free2DMatrix(Y_matrix,line);
    free2DMatrix(res_matrix, DEGREE);
    free2DMatrix(X_matrix_T, DEGREE);
    free2DMatrix(Temp1_matrix, DEGREE);
    free2DMatrix(Temp2_matrix, DEGREE);
}

//预测函数
double toPredict(double x){
    double res=0;
    for (int i = 0; i < DEGREE; ++i) {
        //printf("%lf\n",a[i]);
        res=res+a[i]*pow(x,i+1);
    }
    return res;
}

//生成随机的10个二元对
BinaryPair *initrandomBainay(){
    BinaryPair *head=NULL;
    BinaryPair *p0=NULL;
    BinaryPair *p=NULL;
    int count=1;
    while (count<=DEGREE){
        p=(BinaryPair*) malloc(sizeof (BinaryPair));
        p->next=NULL;
        p->prv=NULL;
        p->x=count;
        //srand(time(NULL));
        p->y=((double)rand())/100;
        if (count==1){
            head=p;
        } else{
            p->prv=p0;
            p0->next=p;
        }
        ++count;
        p0=p;
    }
    return head;
}